Human dermal fibroblast subpopulations and epithelial mesenchymal transition signals in hidradenitis suppurativa tunnels are normalized by spleen tyrosine kinase antagonism in vivo

Hidradenitis Suppurativa is a chronic inflammatory disease of which the pathogenesis is incompletely understood. Dermal fibroblasts have been previously identified as a major source of inflammatory cytokines, however information pertaining to the characteristics of subpopulations of fibroblasts in HS remains unexplored. Using in silico-deconvolution of whole-tissue RNAseq, Nanostring gene expression panels and confirmatory immunohistochemistry we identified fibroblast subpopulations in HS tissue and their relationship to disease severity and lesion morphology. Gene signatures of SFRP2+ fibroblast subsets were increased in lesional tissue, with gene signatures of SFRP1+ fibroblast subsets decreased. SFRP2+ and CXCL12+ fibroblast numbers, measured by IHC, were increased in HS tissue, with greater numbers associated with epithelialized tunnels and Hurley Stage 3 disease. Pro-inflammatory CXCL12+ fibroblasts were also increased, with reductions in SFRP1+ fibroblasts compared to healthy controls. Evidence of Epithelial Mesenchymal Transition was seen via altered gene expression of SNAI2 and altered protein expression of ZEB1, TWIST1, Snail/Slug, E-Cadherin and N-Cadherin in HS lesional tissue. The greatest dysregulation of EMT associated proteins was seen in biopsies containing epithelialized tunnels. The use of the oral Spleen tyrosine Kinase inhibitor Fostamatinib significantly reduced expression of genes associated with chronic inflammation, fibroblast proliferation and migration suggesting a potential role for targeting fibroblast activity in HS.


Introduction
Hidradenitis Suppurativa (HS) is a chronic inflammatory skin disease with an estimated prevalence of 1% globally (Sabat et al., 2020;Wolk et al., 2020) [1,2].Lesions include nodules and abscesses as well as epithelialized dermal tunnels and significant hypertrophic scarring and contractures in flexural areas (Frew et al., 2021a) [3].The molecular pathogenesis of HS is poorly understood, but dermal fibroblasts are known to be major contributors of proinflammatory cytokines including IL-1, CCL20, CXCL1, CXCL8, IL36B and SERPINB4 (Witte-Handel et al., 2018; Zouboulis et al., 2020) [4,5].Fibroblasts are important drivers of wound healing and epidermal repair (Artuc et al., 2002) [6].The mechanisms of wound healing involve activation and migration of fibroblasts as well as epithelial-mesenchymal transition (EMT) like mechanisms which are proposed to be dysregulated in HS (Sen and Roy., 2021; Artuc et al., 2002) [6,7].It has been hypothesized that fibroblasts may play a role in the development of epithelialized tunnels through keratinocyte-fibroblast interactions and EMT (Frew et al., 2019a) [8].The evidence for activation of EMT in lesions of HS is limited and further exploration of this hypothesis may lead to novel therapeutic targets against the progression of disease (Nelson et al., 2019;Frew et al., 2021b) [9,10] including the development of tunnels, hypertrophic scarring and flexural contractures.
The therapeutic modulation of fibroblasts is an emerging area of research in fibroinflammatory diseases with existing therapies including Janus Kinase (JAK) inhibitors and Spleen Tyrosine Kinase (SYK) inhibitors demonstrating evidence of clinically relevant fibroblast modulation (Liu et al., 2019) [16].These drugs have the potential to be repurposed for investigation into their benefit in the setting of HS.SYK inhibition has the additional benefit of potentially addressing B cells, plasma cells and monocytes which are significantly dysregulated in HS (Gudjonsson et al., 2020, Musilova et al., 2020, Byrd et al., 2019;Frew et al 2020) [17][18][19][20].
In this study we aimed to examine the presence of fibroblast subpopulations in lesional tissue (Frew 2019b) [21] of HS compared to non lesional tissue and healthy controls.Additionally, the effect of SYK antagonism with Fostamatinib upon gene expression profiles in lesional tissue was assessed.

Immunohistochemistry validates the presence of discrete fibroblast subpopulation locations in hidradenitis suppurativa tissue
Transcriptomic findings illustrating differential gene expression of fibroblast sub-populations was validated using immunohistochemistry with primary antibodies against SFRP1, SFRP2 and CXCL12 proteins (primary antibodies used listed in S3 Table) in lesional and non lesional HS tissue (Fig 3a -3c).Significant differences in the number of SFRP1+, SFRP2+ and CXCL12 + cells (measured using semiquantitative IHC staining) was seen between normal and HS lesional tissue (Fig 3d -3f).Stratification by Hurley stage illustrated a significant reduction in SFRP1+ cells and an elevation in the number of SFRP2+ and CXCL12+ cells in Hurley Stage 3 patients compared with Hurley stage 2 patients (Fig 3d -3f).Stratification by the presence of epithelialized tunnels (Navrazhina et al., 2021) [23] illustrated no significant difference between the number of SFRP1+ cells, but a decrease in the number of SFRP2+ cells and CXCL12+ cells in the absence of epithelialized tunnels compared to the presence of epithelialized tunnels (Fig 3d -3f).Differences in intensity of N-Cadherin staining were observed between the overlying epidermal epithelium when compared with the epithelial lining of tunnels (Fig 4h -4k).Significant increases in the number of N-Cadherin positive cells were seen in the overlying epidermal epithelium than epithelialised tunnels.The semiquantitative IHC staining in epithelialized tunnels was comparable to that of healthy controls (Fig 4j).Conversely, E-Cadherin protein expression is negligible in overlying epidermis but increased in tunnel epithelium compared to healthy controls, however this difference did not reach statistical significance (Fig 4k).

Spleen tyrosine kinase inhibition with fostamatinib demonstrates significant reduction in fibrosis, inflammatory and epithelial-mesenchymal transition associated genes and pathways in hidradenitis suppurativa after 4 weeks of therapy
Spleen Tyrosine Kinase (SYK) inhibition with Fostamatinib at a dose of 100mg twice daily significantly modulates gene expression as measured by the Human Fibrosis V2.0 Nanostring

Discussion
This study presents novel gene expression and immunohistochemical data pertaining to the presence of fibroblast subpopulations and evidence of EMT in lesional tissue of HS.Fibroblast subpopulations in lesional tissue are significantly altered compared to non lesional tissue and healthy controls.Expansion of the SFRP2+ and CXCL12+ populations were identified with no significant alteration to the SFRP1+ population.Additionally, the severity of disease and presence of epithelialized tunnels is associated with these alterations to fibroblast subpopulations.Evidence of EMT was seen in HS lesional tissue with greater levels associated with epithelialised tunnels and Hurley stage 3 disease.SYK inhibition with fostamatinib for as little as 4 weeks resulted in reduced expression of fibroblast and EMT associated genes as measured by Nanostring gene expression assay.Pathway analysis of Nanostring data demonstrated enrichment of downregulated pathways pertaining to fibroblast proliferation and connective tissue growth factor production.Previously published functional annotation of fibroblast subpopulations in normal human skin identifies SFRP2+ (Type A) cells as functionally associated with dermal cell and extracellular matrix homeostasis as well as inflammatory cell retention.Additionally, SFRP2+ cells are involved in the development of tertiary lymphoid structures (TLO) (Tabib et al., 2018) [14].The significant expansion of SFRP2+ fibroblast subpopulations in HS lesional tissue would explain the development of TLOs (Wolk et al., 2020;Frew et al., 2020) [2,18] in HS tissue.Additionally, gross histological findings in HS include gross expansion of the dermal layer, fibrosis and retention of inflammatory cells, all associated with function of SFRP2+ fibroblasts.
CXCL12+ (Type B) cells contribute to immune surveillance and are associated with Th2 immune axis activation.This subset is particularly expanded in lesional tissue of atopic dermatitis (He et al., 2019) [24].SFRP1+ cells are described as being functionally responsible for ECM remodeling and include dermal papilla cells and fibroblasts of the dermal-hypodermal junction.The reduction in SFRP1+ cells in advanced Hurley stage 3 disease may be reflective of destruction of follicular units and decreases in the number of dermal papilla fibroblasts in tissue.Additionally, CXCL12+ cells were increased in severe disease and the presence of epithelialized tunnels suggesting an association with increased levels of inflammation (Navrazhina et al.,2021) [23].Witte Handel et al (Witte Handel et al., 2018) [4] has previously demonstrated fibroblasts as a major source of IL-1 and downstream inflammasome activation in HS, and our data supports the hypothesis that SFRP2+ and CXCL12+ fibroblasts are the subpopulations associated with this inflammatory signature.Further investigations with single cell RNA sequencing would be required to confirm this, as tissue culture results in preferential expansion of SFRP2+ cells (Tabib et al., 2021) [15] which may lead to poor isolation of other fibroblast subsets.
Evidence of EMT has been reported in HS lesional tissue by Nelson and colleagues with reductions in E-cadherin in the epidermis (Nelson et al., 2019) [10].This has been interpreted as a sign of sensitivity to mechanical stress (Danby et al., 2013) [25], however our investigations in the epithelium of tunnels suggest that greater EMT signals are associated with epithelialised tunnels and may be indicative of a mechanism of tunnel development rather than a predisposing factor of disease initiation (Frew et al., 2019a) [8].This is supported by the fact that nonlesional tissue did not demonstrate any significant differences in EMT protein expression compared to normal skin (Fig 4a -4d), and that the EMT expression as measured by IHC was significantly elevated in biopsies with tunnels compared to those without (Fig 4e -4g).
EMT can be induced in epithelial cells by inflammation and extracellular matrix alterations in both oncological and inflammatory settings (Suarez-Carmona., 2017; Wang et al., 2018) [26,27].Additionally, the process of EMT stimulated the production of inflammatory cytokines including IL-1 (Suarez-Carmona., 2017) [26].These results clearly demonstrate EMT signals in tissue of HS, however the precise causative cascade leading to the development of EMT in HS remains unclear.The use of Fostamatinib partially ameliorates the gene expression signature of EMT in HS tissue, presumably indirectly through its action upon B cells, plasma cells and monocyte/macrophages (Frew., 2020; Wang et al., 2021) [18,28].Wang demonstrates direct alteration of fibroblast function with the use of fostamatinib in vitro (Wang et al., 2021) [28] and therefore fostamatinib may have potential as an adjuvant treatment alongside other targeted therapies such as Janus Kinase Inhibitors or Monoclonal antibodies (Frew et al., 2021b) [9], particularly in the setting of epithelialised tunnels and/or significant fibrosis.
The use of whole genome RNAseq has several disadvantages when examining discrete gene expression pertaining to cellular subsets (Whitley et al., 2016) [29].Biases during cDNA library construction and sequence alignment, as well as read depth can alter transcript quantitation and reduce reproducibility between datasets.Single cell RNA sequencing is susceptible to drop out and significant noise bias due to the low amount of RNA analyzed (Wu et al., 2018) [30].
Nanostring gene expression arrays addresses some of these disadvantages including no need for amplification, ability to identify low volume and degraded RNA and minimal background signal compared to whole tissue RNAseq.When RNAseq and Nanostring gene expression analyses have been compared in Psoriasis Vulgaris, increased accuracy for low level expressed transcripts was identified in Nanostring compared with whole genome RNAseq (Krueger et al., 2019) [31].
In-silico deconvolution provides a computational alternative to scRNAseq which, when utilized with a high-quality reference frame can provide insights into cellular proportions in whole RNAseq data whilst avoiding the biases present in scRNAseq.The quality of the deconvolution data is predicated on the baseline RNAseq data.Whilst whole tissue RNAseq is most appropriate for biomarker discovery, for identification of a priori differential gene expression, the combination of in silico-deconvolution of whole tissue RNAseq and Nanostring gene expression (confirmed with immunohistochemistry) provides a robust methodology for investigation of pre-defined signatures of fibroblast subpopulations.
HS is in great need of novel therapies; however, our understanding of the disease pathogenesis and specific mechanistic action of novel therapies is incomplete.SYK is an important mediator of downstream signaling in monocytes, neutrophils and B cells but can also contribute to the reprogramming of fibroblasts (Liu et al 2019) [16].The positive clinical results in our recent open-label clinical trial of fostamatinib in HS (Jepsen et al., 2023) [32] suggest that fostamatinib has a clinically significant impact on disease severity in HS.The analysis of gene expression data from the first 4 weeks of therapy suggests that a rapid alteration in chronic inflammatory and pro-fibrotic pathways is associated with this therapy and is associated with clinical response.Future analysis of this clinical trial data will enable deeper investigation into the precise mechanisms of SYK inhibition in HS at higher doses and further timepoints.

Patients and tissue collection
This study was approved by the human research ethics committee of Sydney South West Area Health Service (Sydney Australia).A total of 20 participants underwent lesional and non lesional biopsies (6mm in diameter) based upon previously published definitions (Frew 2019b) [21].All biopsies were immediately stored in RNA Later and frozen at -80C until processing.Standard wound care was undertaken for each of the biopsy sites.All patients were not on any active therapy at the time of biopsies.The clinical characteristics of participants are included in S1 Table.

Skin biopsy mRNA quantification
Total RNA was extracted using the Qiagen RNeasy kit as per standard protocol (Qiagen, Valencia, CA, USA).RNA was sequenced using the Novaseq S4 (Illumina, San Diego, CA, USA) Differential expression was estimated with DESeq2 and VST to log 2 transform the normalized counts.Differentially expressed genes were defined as fold change (FCH) of � 1�5 and false discovery rate � 0�05.

Pathway analysis
Enrichment analyses were completed using Protein ANalysis THrough Evolutionary Relationships (PANTHER) classification system (www.pantherdb.org).Analyses were performed using whole tissue RNAseq data (

CIBERSORTx
In silico deconvolution was undertaken using the CIBERSORTx web application ((https:// cibersortx.stanford.edu/index.php) in order to estimate the abundance of fibroblast subpopulations in baseline lesional and non-lesional tissue using whole tissue RNAseq.The utilized fibroblast reference matrix was based upon previously published human skin scRNAseq data (GSE128169) with stratification into three major subtypes as defined by Ascension et al (https://doi.org/10.1016/j.jid.2020.11.028).

Nanostring analysis
Extracted RNA was analysed using the NCounter system (Nanostring) using the Human Fibrosis V2.0 gene panel.Raw data was processed using the Nanostring Nsolver (version 4.0.70)analysis software using quality control and normalization procedures derived from the NormqPCR R package as previously described (Bhattacharya et al., 2021) [33].Differentially expressed genes (DEGs) were defined as >1.5 Log2Fold change with a false discovery rate<0.05.

Immunohistochemistry
IHC staining was undertaken using primary antibodies listed in S3 Table .Stained slides underwent semiquantitative analysis using standardized, previously published methods (Crowe et al., 2019) [34].Quantitation of positive staining cells was performed by mm 2 of epidermis or dermis and measured using ImageJ v1.42 software.Differences between groups were analyzed using the Wilcoxon rank sum test for two variables and/or omnibus test using one-way ANOVA with adjustment for multiple comparisons was made using the Benjamini Hochberg procedure.

Fostamatinib administration
Fostamatinib was administered as part of the Phase 2 clinical trial in moderate to severe HS (NCT0504698) (Clinicaltrials.gov)[35].100mg BID was administered for 4 weeks with 150mg BID thereafter until week 12 of therapy.Specific details regarding drug administration, monitoring and clinical outcomes are published separately.

Statistical analysis
RNA sequencing statistical analysis was performed in R language (www.r-project.org;R Foundation, Vienna, Austria) using publicly available packages (www.bioconductor.org).Means of each group and differences between groups were estimated using least-square means.Hypothesis testing was undertaken using the general framework for linear mixed-effect models of the limma package.tissue was considered a fixed factors while random effect related to subjects was included.Unsupervised clustering was performed with the Euclidean distance and average agglomeration.RNAseq data p-values were corrected for multiplicity using Benjamini-Hochberg procedure.

Fig 2 .
Fig 2. In-Silico Deconvolution of whole tissue RNAseq demonstrates differential expansion of fibroblast subpopulations in lesional vs non lesional tissue of Hidradenitis Suppurativa.(A): Cell Proportion data as calculated by CIbersortX in-silico deconvolution demonstrating alterations in cell populations between lesional and non lesional tissue.(B) Isolation of only fibroblast subsets indicates an expansion of SFRP2+ fibroblast subpopulation gene signatures.(C): Individual sample data identifies a statistically significant expansion of SFRP2+ fibroblast subsets in lesional tissue, decreased SFRP1+ populations and no significant change in CXCL12+ populations in whole tissue RNAseq with in-silico deconvolution.https://doi.org/10.1371/journal.pone.0282763.g002

Fig 3 .
Fig 3. (A-C) Fibroblast subpopulations in HS differ based upon Hurley stage and presence of epithelialised tunnels.SFRP1 protein (D) is diffusely positive in the reticular dermis in Hurley stage 2 lesions but reduced in Hurley 3. No difference was seen in the presence of tunnels.SFRP2 protein (E) illustrates increased expression in Hurley stage 3 disease and in the presence of epithelialised tunnels.CXCL12 protein expression (F) is increased in Hurley stage 3 disease and in the presence of epithelialized tunnels.Length of marker equals 1mm.https://doi.org/10.1371/journal.pone.0282763.g003 Immunohistochemical staining for signs of epithelial mesenchymal transition (EMT) identified increased staining of TWIST1 (Fig 4a), ZEB1 (Fig 4b), Snail/Slug (Fig 4c) and N-Cadherin (Fig 4d) proteins in HS lesional tissue compared to site-matched healthy controls.The tips of elongated epithelial rete ridges in pseudo-psoriasiform hyperplasia demonstrate discrete clustering of positive staining to EMT-associated proteins including TWIST1 and ZEB1 (Fig 4a and 4b).Stratification of HS lesional tissue by the presence of epithelialized tunnels indicated a significant increase in the number of SLUG/SNAIL, ZEB1 and TWIST1 positive cells when compared with HS lesions without tunnels or healthy controls.(Fig 4e-4g).

Fig 4 .
Fig 4. Epithelial Mesenchymal Transition Markers are significantly overexpressed in HS lesional tissue compared to site matched healthy controls (A-D).Protein expression of TWIST1 (E), ZEB1 (F) and SLUG/SNAIL (G) are significantly increased in HS lesional tissue with highest expression in lesions with epithelialized tunnels.Protein expression of N Cadherin (H) was significantly elevated in the overlying epidermis of lesions with tunnels but normal in the epithelial lining of the tunnel itself (J).E-Cadherin expression (I) showed no significant difference between sites (K).Length of marker equals 1mm.https://doi.org/10.1371/journal.pone.0282763.g004

Fig 5 .
Fig 5. Therapy with 4 weeks of the oral spleen tyrosine kinase antagonist Fostamatinib downregulates gene expression relating to fibroblast activity and epithelial-mesenchymal transition as measured by Nanostring gene expression panel (A,B).Pathway analysis identified significant downregulation and enrichment of chronic inflammatory response, connective tissue growth factor and transforming growth factor beta production (C).https://doi.org/10.1371/journal.pone.0282763.g005 Fig 1) as well as for Nanostring Gene Expression Panels.Nanostring gene expression panels included Baseline (prior to fostamatinib administration) and Week 4 of fostamatinib therapy lesional tissue samples (Fig 5).